CNN evaluation
#rm(list = ls())
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(future)
library(ggthemes)In this notebook we will test the performance of varKode to distinguish species of Malpighiaceae and figure out the best parameters for training a dataset.
Kmer size and amount of data
To start, we produced images from different numbers of kmers. We can suppose that shorter kmers will offer lower resolution to resolve species, but they will also create smaller files that require less computation. Here we will test whether images based on longer kmers result in higher accuracy. As an example, here are images produced from 200Mb for the same sample, but different kmer sizes (5-9):
knitr::include_graphics(paste0('images_',5:9,'/S_bannisterioides+S-91_00200000K.png'))
We also used different amounts of data to produce images, since we want
to figure out the lowest amount needed to distinguish species. With less
data, figures get more noisy since chance plays a bigger role in the
observed kmer frequencies. This should be more severe for larger kmer
sizes, since each kmer will be more unique in the genome.
For example, images for 5-mer for the same sample as above, for 500Kb and 200Mb:
knitr::include_graphics(paste0('images_6/S_bannisterioides+S-91_00',c('000500','200000'),'K.png'))
The same, but for 8-mers:
knitr::include_graphics(paste0('images_8/S_bannisterioides+S-91_00',c('000500','200000'),'K.png'))
Now that we understand the differences between images, let’s understand the effect in accuracy. We previously trained CNN models to recognize images for a combination of kmer sizes and amount of data, with 10 replicates for each combination. In each replicate, we kept 3 randomly chosen samples per species as a validation set and checked the accuracy of the trained model in guessing the species of these samples, for different amounts of data used for the validation sample. What we want is to find:
1 - The lowest kmer size to produce high accuracy
2 - The lowest amount of data needed
3 - Whether the amount of data used for training and for querying must be similar.
The results of these simulations were saved as a csv table, let’s load it (ignoring the first, index column):
df = read_csv('kmerSize_VS_bp.csv')[-1]## New names:
## Rows: 4500 Columns: 11
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (3): bp_training, samples_training, samples_valid dbl (8): ...1, kmer_size,
## replicate, bp_valid, n_samp_training, n_samp_valid...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`dfNow let’s make sure bp_training and bp_valid are treated as ordered factors for nice plotting:
not_all = as.character(sort(as.integer(unique(df$bp_training[!str_detect(df$bp_training,'\\|',)]))/1e6))
ordered_levels = c(not_all,'all')
df = df %>%
mutate(bp_training = as.character(as.integer(bp_training)/1e6) ) %>%
mutate(bp_training = replace_na(bp_training, 'all')) %>%
mutate(bp_training = factor(bp_training,
levels = ordered_levels,
ordered = TRUE),
bp_valid = factor(as.character(as.integer(bp_valid)/1e6),
levels=ordered_levels,
ordered = TRUE),
kmer_size = factor(as.character(kmer_size),
levels = as.character(sort(unique(kmer_size))),
ordered = TRUE
)
)## Warning in mask$eval_all_mutate(quo): NAs introduced by coerciondfNow we can plot:
kmer_labeller = as_labeller(function(value){
return(paste0('kmer length:',value))
})
ggplot(df) +
geom_jitter(aes(x = bp_training, y = bp_valid, color = valid_acc)) +
scale_color_viridis_c('Validation\naccuracy', option = 'inferno', limits = c(0,1)) +
facet_grid(~kmer_size, labeller = kmer_labeller) +
coord_equal() +
xlab('Data in training images (Mb)') +
ylab('Data in validation images (Mb)')
Now a version with averaged accuracy
p = df %>%
group_by(kmer_size,bp_training,bp_valid) %>%
summarize(valid_acc = mean(valid_acc)) %>%
ggplot(aes(x = bp_training, y = bp_valid, fill = valid_acc)) +
geom_raster() +
#geom_text(aes(label=sprintf(100*valid_acc,fmt='%2.0f')),size=4.5*5/14) +
scale_fill_viridis_c('Average\nvalidation\naccuracy', option = 'magma', limits = c(0,1),labels=scales::percent) +
facet_grid(~kmer_size, labeller = kmer_labeller) +
coord_equal() +
xlab('Data in training images (Mb)') +
ylab('Data in validation images (Mb)') +
theme_few(base_size = 6)## `summarise()` has grouped output by 'kmer_size', 'bp_training'. You can
## override using the `.groups` argument.p
dir.create('paper_images')## Warning in dir.create("paper_images"): 'paper_images' already existsggsave(filename = 'kmerlen_vs_accuracy.png',plot =p,device='png',path = 'paper_images',width = 22,height = 5,units = 'cm',dpi = 2400)means = df %>%
filter(bp_training %in% c('0.5','1','200','all')) %>%
#filter(bp_valid %in% c('50','100','200','all')) %>%
filter(bp_valid %in% c('2','5','10','20','50','100','200')) %>%
group_by(bp_training,kmer_size) %>%
summarise(Int=median(valid_acc))## `summarise()` has grouped output by 'bp_training'. You can override using the
## `.groups` argument.df %>%
filter(bp_training %in% c('0.5','1','200','all')) %>%
filter(bp_valid %in% c('2','5','10','20','50','100','200')) %>%
#filter(bp_valid %in% c('50','100','200','all')) %>%
ggplot(aes(x=valid_acc)) +
geom_histogram(aes(x=valid_acc)) +
facet_grid(kmer_size~bp_training) +
geom_vline(data = means, aes(xintercept = Int))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
So it seems that the smallest kmer sizes never result in very high accuracy, and the largest kmer sizes result in high accuracy for higher amounts of data, but lower accuracy for lower amounts. It seems that a kmer size of 7 is a good balance, and that training using images of different sizes helps in being more robust to the amount of data used to produce validation images.
As little as 1Mb produces moderately accurate results for kmer size 7 or below.
Can we quantify what is different about images produced with different data amounts?
images = c(list.files(path='images_5',pattern='.png', recursive = T, full.names = T),
list.files(path='images_6',pattern='.png', recursive = T, full.names = T),
list.files(path='images_7',pattern='.png', recursive = T, full.names = T),
list.files(path='images_8',pattern='.png', recursive = T, full.names = T),
list.files(path='images_9',pattern='.png', recursive = T, full.names = T))
nkmers = function(k){ #from https://bioinfologics.github.io/post/2018/09/17/k-mer-counting-part-i-introduction/
(4^k + (1 - k%%2) * 4^(k/2))/2
}
get_sd = function(path){
k = as.integer(gsub('.+_([0-9])/.+','\\1', path))
taxon = gsub('.+/(.+)\\+.+','\\1', path)
sample = gsub('.+\\+(S-[0-9]+)_.+','\\1', path)
Mbp = as.integer(gsub('.+_([0-9]{8})K.+','\\1', path)) / 1000
x = sort(png::readPNG(path))
x = x[(length(x)-nkmers(k)+1):length(x)]
sd_counts = sd(table(x))
data.frame(k = k, taxon = taxon, sample = sample, Mbp = Mbp, sd_counts=sd_counts)
}
plan(multisession(workers = 4))
df = furrr::future_map_dfr(images,get_sd)
plan(sequential)
dfggplot(df) +
geom_line(aes(x=Mbp, y=sd_counts,color=sample)) +
facet_wrap(as.factor(k)~.,scales = 'free') +
scale_color_discrete(guide='none') +
scale_y_log10() +
scale_x_log10()
Training parameters
Now we will check the results of using different training parameters: - model pretraining - augmentation (CutMix or MixUp) - Label Smoothing - model architecture - lighting transformations
Let’s read the data and prepare for plotting:
df = read_csv('training_params.csv')[-1] %>%
mutate(bp_valid = factor(as.character(as.integer(bp_valid)/1e6),
levels = sort(unique(bp_valid/1e6)),
ordered = TRUE),
augmentation = ifelse(str_detect(callback,'CutMix'),'CutMix',
ifelse(str_detect(callback,'MixUp'),'MixUp',
'None')
),
augmentation = factor(augmentation, levels = c('None','MixUp','CutMix'),ordered = F),
aug = str_replace(augmentation,'None',''),
lablsmth= ifelse(label_smoothing,
'label Smoothing',
''),
pretr = ifelse(pretrained,
'pretrained',
''
),
transformations = ifelse(trans,
'with_transforms',
''
),
parameters = paste(arch,pretr,lablsmth,aug,transformations,sep=',') %>%
str_replace_all(',{2,}',',') %>%
str_remove_all('^,|,$') %>%
str_replace_all('^$','None') %>%
fct_reorder(valid_acc, mean)
)## New names:
## Rows: 30240 Columns: 16
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (5): bp_training, samples_training, samples_valid, callback, arch dbl (8):
## ...1, kmer_size, replicate, bp_valid, n_samp_training, n_samp_valid... lgl (3):
## label_smoothing, pretrained, trans
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`df Now we can plot the effect of parameters. There are clearly some models that do much better than others:
ggplot(df, aes(x = parameters, y = valid_acc)) +
#geom_boxplot() +
#geom_violin(adjust=1.5) +
geom_jitter(aes(color = bp_valid),height = 0.005) +
scale_color_viridis_d(option='turbo',begin = 0.1, end=0.9) +
#facet_wrap(~bp_valid) +
theme(axis.text.x = element_text(hjust = 1, angle = 45))
What are the 20 most accurate models?
ggplot(filter(df, parameters %in% tail(levels(df$parameters),20)), aes(x = parameters, y = valid_acc)) +
#geom_boxplot() +
#geom_violin(adjust=1.5) +
geom_jitter(aes(color = bp_valid),height = 0.005) +
scale_color_viridis_d(option='turbo',begin = 0.1, end=0.9) +
#facet_wrap(~bp_valid) +
theme(axis.text.x = element_text(hjust = 1, angle = 45))